cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar
import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
from scipy.stats import boxcox
clima = gpd.read_file('climas_hemispheres')
clima.columns
area = clima['SUM_Area']*1e-6
data_columns = clima.columns[3:3+17*12]
data_columns
inter_cols = {}
for yr in range(17):
inter_cols[yr] = data_columns[12*yr:12*yr +12]
intra_cols = {}
for mo in range(12):
intra_cols[mo] = [data_columns[mo + yr *12 ] for yr in range(17)]
for yr in range(17):
inter = clima[ inter_cols[yr]].sum(axis=1)/area
label = 'inter_' + str(yr +2001)
clima= clima.assign(label=inter)
clima = clima.rename(columns={'label': label})
inter_columns = clima.columns[-17:]
inter_columns
inter_max = clima[inter_columns].max().max()
for mo in range(12):
intra = clima[ intra_cols[mo]].sum(axis=1)/(17*area)
label = 'intra_{:02}'.format(mo)
clima = clima.assign(label=intra)
clima = clima.rename(columns={'label': label})
intra_columns = clima.columns[-12:]
intra_columns
intra_max = clima[intra_columns].max().max()
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100,cmap='OrRd'):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap=cmap,vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=8)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
# title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in ' +mo_name + ' (2001-2017)'
lambdas = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column, lbd = boxcox(clima[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
lambdas
maxs = []
mins = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column = boxcox(clima[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
for yr in range(17):
column = 'inter_{}'.format(yr +2001)
test_column = (boxcox(clima[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
clima = clima.assign(test_column =test_column)
yr_name = str(yr +2001)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +yr_name
name = 'inter_climas_ns_index_{}'.format(yr +2001)
draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name)
lambdas = []
for mo in range(12):
column= 'intra_{:02}'.format(mo)
test_column, lbd = boxcox(clima[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
lambdas
maxs = []
mins = []
for mo in range(12):
column= 'intra_{:02}'.format(mo)
test_column = boxcox(clima[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
for mo in range(12):
column = 'intra_{:02}'.format(mo)
test_column = (boxcox(clima[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
clima = clima.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +mo_name + ' (2001-2017)'
name = 'intra_climas_index_{:02}'.format(mo)
draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name)
total = clima[data_columns].sum(axis=1)/17
total = total/area
clima= clima.assign(total=total)
column = 'total'
test_column, lbd = boxcox(clima[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
clima = clima.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere (2001-2017)'
name = 'total_climas_nsparts_index'
draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name, dpi=600)
inter_medias = clima[inter_columns].mean(axis=1)
clima[inter_columns].shape
anomalias = np.array(clima[inter_columns]) - inter_medias.values.reshape(49,1)
bigger = abs(anomalias.max()) if abs(anomalias.max()) >= abs(anomalias.min()) else abs(anomalias.min())
vmax = bigger; vmin = -1 * bigger
for yr in range(17):
test_column = anomalias[:,yr]
clima = clima.assign(test_column =test_column)
yr_name = str(yr +2001)
title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +yr_name
name = 'inter_anomalies_climas_ns_{}'.format(yr +2001)
draw_map(clima,'test_column',title=title,name=name,cmap='seismic',vmin=vmin,vmax=vmax)
intra_medias = clima[intra_columns].mean(axis=1)
anomalias = np.array(clima[intra_columns]) - intra_medias.values.reshape(49,1)
bigger = abs(anomalias.max()) if abs(anomalias.max()) >= abs(anomalias.min()) else abs(anomalias.min())
vmax = bigger; vmin = -1 * bigger
for mo in range(12):
test_column = anomalias[:,mo]
clima= clima.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +mo_name + ' (2001-2017)'
name = 'intra_anomalies_climate_ns_{:02}'.format(mo)
draw_map(clima,'test_column',title=title,name=name,cmap='seismic',vmin=vmin,vmax=vmax, dpi =600)